Canek integration analysis of TIL-NK (WT/KO)

library(tidyverse)
library(Seurat)
library(Canek)
library(scmisc)
set.seed(1)
theme_set(cowplot::theme_cowplot())

Load data

y <- read_rds("~/NK_integration_process_livenk_sampling.rds")

DefaultAssay(y) <- "RNA"

DimPlot(y, split.by="background", label = TRUE) + NoLegend()+
   ggtitle("Canek integration of WT and KO NK cells (TIL)")+
  theme(plot.title = element_text(size=13, face = "bold.italic"))

VlnPlot(y, c("nCount_RNA", "nFeature_RNA", "mito.percentage"), pt.size = .02)

k <- y[,!y$seurat_clusters == "11"]
DimPlot(k, label=TRUE)

set.seed(1)

k <- FindVariableFeatures(k, nfeatures = 5000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

VariableFeaturePlot(k)

k <- ScaleData(k)

Centering and scaling data matrix

  |                                                                                                    
  |                                                                                              |   0%
  |                                                                                                    
  |===================                                                                           |  20%
  |                                                                                                    
  |======================================                                                        |  40%
  |                                                                                                    
  |========================================================                                      |  60%
  |                                                                                                    
  |===========================================================================                   |  80%
  |                                                                                                    
  |==============================================================================================| 100%

k <- RunPCA(k)

PC_ 1 
Positive:  Ccl5, Dusp1, Cd7, Gm42418, Klf2, Zfp36l1, Zfp36, Fos, Emb, Pmepa1 
       Il18r1, Gm26740, Klf3, Nr4a2, Klrb1b, Sult2b1, Irf7, Nfkbia, Jun, Ltb 
       Sat1, Nr4a1, Lmna, Cd28, Tcf7, Klra8, Gpc1, Klra4, Cited2, Dhrs3 
Negative:  Pclaf, Dut, Cks1b, Spc24, Srm, Ran, Ranbp1, Birc5, Tk1, Phgdh 
       C1qbp, Nme1, Nme2, Tyms, Txn1, Mif, Ppia, Rrm2, Rps27l, Mcm3 
       Nhp2, Rps17, Ldha, Nap1l1, Mcm5, Ccna2, Stmn1, Clspn, Mki67, Snrpd1 
PC_ 2 
Positive:  Rpl36, Rps26, Cd52, Rpl37, Rpl38, Rpl35, Rpl41, Rpl39, Rpl32, Rps28 
       Rps8, Rpl27, Nme1, Sec61g, Serp1, Rpl37a, Rps29, Atp5e, Nme2, Tmsb4x 
       Ehd1, Sec61b, Rplp2, Srm, Ctla4, Rpl28, Ppp1r14b, Rpl10a, Rplp0, Hilpda 
Negative:  Ccna2, Pclaf, Mki67, Rrm2, Kif11, Esco2, Birc5, Aurkb, Hist1h1b, Nusap1 
       Cdca8, Tpx2, Hist1h2ae, Top2a, Spc24, Fbxo5, Tk1, Stmn1, Clspn, Cdk1 
       Cdca3, Asf1b, Cenpf, Cenpe, Kif4, Rad51, Neil3, Knl1, Kif22, Shcbp1 
PC_ 3 
Positive:  Ptma, Eif5a, Xcl1, Fbl, Npm1, Srm, Mybbp1a, Hspa8, Ccnd2, Emb 
       Eif4a1, Cd82, Hsp90ab1, Nr4a2, Ncl, Impdh2, Bzw2, Srsf2, Hspa9, Nolc1 
       Mat2a, Gnl3, Nop58, Snhg1, Hnrnpab, Ppp1r14b, C1qbp, Pa2g4, Ranbp1, Rpl14 
Negative:  Klrg1, S100a4, Cxcr6, Epas1, Irak3, S100a6, Glrx, Ctla4, Adgrg1, Ly6a 
       Tmsb4x, Cd3g, Itga1, Gzmc, S100a11, Filip1, Gng2, Entpd1, Serpinb6a, Hilpda 
       AA467197, Sema6d, Plac8, Vim, Tmem154, Fgl2, Fcgr3, Ebi3, Tnfrsf13b, Ehd1 
PC_ 4 
Positive:  Rpl39, Rpl37, Rpl32, Rpl37a, Rps20, Rps28, Rpl36, Rpl41, Rps26, Rplp0 
       Rps29, Rpl38, Rpl36a, Rpl10a, Rps12, Rps19, Rpl35, Rplp2, Rps8, Rpl28 
       Rps18, Rpl26, Rpl27, Emb, Rpl12, Rpl15, Tcf7, H2-Oa, Rps2, Kit 
Negative:  Gzma, Gzmb, Ccl5, Prf1, Actb, Ccl4, Emp3, Cyba, Itgb2, Irf8 
       Pfn1, Itgam, Tnfrsf9, Klrg1, Pim1, Lgals1, Icam1, Lcp1, Epas1, Adgrg1 
       Ctsd, Ccnd2, Entpd1, Plek, Bcl2a1b, Actg1, Calm1, Gzmc, Ccl3, Cma1 
PC_ 5 
Positive:  Gapdh, Rpl15, Rplp0, Rps2, Nrgn, Rps18, Lilrb4a, Ccr2, Pglyrp1, Rpl10a 
       Ly6e, Rps20, Emb, Ier5l, Rpl7a, Xist, Rps12, Nr4a2, Gpi1, Gm43305 
       Thy1, Rbm3, Cd82, Rpl14, Rpl12, Rps6, Atp5g2, Anxa1, Lilr4b, Rgs2 
Negative:  Gm42418, Cd3e, Pou2f2, Crip1, Zc3h12a, Nfkbiz, Acp5, Srm, Rpl38, Serp1 
       Rps29, Rpl37a, Cd3d, Rpl36, Eif2s3y, Rilpl2, Rasgrp2, Hspa1b, Nolc1, Rpl37 
       Ccl3, Ifng, S100a6, Hspd1, Hspa1a, Gimap7, Klf2, Cma1, Dus2, Irak3 

ElbowPlot(k,ndims = 50)

k <- RunUMAP(k, dim = 1:20)

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session12:18:24 UMAP embedding parameters a = 0.9922 b = 1.112
12:18:24 Read 18929 rows and found 20 numeric columns
12:18:24 Using Annoy for neighbor search, n_neighbors = 30
12:18:24 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
12:18:26 Writing NN index file to temp file /var/folders/gw/069y9h752hdbn01wcy57ktj00000gn/T//RtmpFi8vRo/filefa2719ceb83
12:18:26 Searching Annoy index using 1 thread, search_k = 3000
12:18:31 Annoy recall = 100%
12:18:32 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
12:18:34 Initializing from normalized Laplacian + noise (using irlba)
12:18:34 Commencing optimization for 200 epochs, with 816206 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
12:18:45 Optimization finished

DimPlot(k, label = TRUE) 

k <- FindNeighbors(k, dims = 1:20)

Computing nearest neighbor graph
Computing SNN

k <- FindClusters(k, algorithm = 4, resolution = 0.6)

Fig.S4D

DimPlot(k, split.by="background", label = TRUE) + NoAxes()+
   ggtitle("Integration analysis of WT and KO NK cells from TIL single cell samples")+
  theme(plot.title = element_text(size=13, face = "bold.italic"))

Fig.S4E

Fig.S4F

VlnPlot(k, features = c("Ifng","Gzmb"), group.by = "background", pt.size = 0.0)

Fig.S4I

chemo<- c("Cxcr1","Cxcr2","Cxcr3","Cxcr4","Cxcr5","Cxcr6",
          "Ccr1","Ccr2","Ccr3","Ccr4","Ccr5","Ccr6","Ccr7","Ccr8","Ccr9","Ccr10",
          "Xcr1","Cx3cr1")

DotPlot(k, features = chemo,group.by = "background") + RotatedAxis() +
    scale_colour_distiller(palette = "RdBu")

Warning: Scaling data with a low number of groups may produce misleading resultsScale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

Fig.4D

VlnPlot(k, "Cxcr6", group.by = "background", pt.size = 0.02)

k<- NormalizeData(k)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

deg <- FindAllMarkers(k,max.cells.per.ident = 500, verbose = FALSE)

Fig.S4H

plot_volcano(deg, n=20, fdr=0.05, groups = "10")

top <- deg%>%filter(p_val_adj<0.01, abs(avg_log2FC)>1)

top_genes <- top_deg(deg)

z <-sample_cells(k, group = "seurat_clusters")

Neither 'n', 'frac' or 'n_max' specified. Sampling 52 cells per group.

Fig.S4G

plot_heatmap(z[top_genes$gene,], column_split= z$seurat_clusters)

sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS 13.3.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] scmisc_0.8.1       Canek_0.2.2        SeuratObject_4.1.3 Seurat_4.3.0.1     lubridate_1.9.2    forcats_1.0.0      stringr_1.5.0     
 [8] dplyr_1.1.3        purrr_1.0.2        readr_2.1.4        tidyr_1.3.0        tibble_3.2.1       ggplot2_3.4.3      tidyverse_2.0.0   

loaded via a namespace (and not attached):
  [1] utf8_1.2.3                  spatstat.explore_3.2-1      reticulate_1.31             tidyselect_1.2.0            htmlwidgets_1.6.2          
  [6] grid_4.1.0                  BiocParallel_1.28.3         Rtsne_0.16                  munsell_0.5.0               codetools_0.2-19           
 [11] ica_1.0-3                   future_1.33.0               miniUI_0.1.1.1              withr_2.5.0                 spatstat.random_3.1-5      
 [16] colorspace_2.1-0            progressr_0.14.0            Biobase_2.54.0              knitr_1.43                  rstudioapi_0.15.0          
 [21] stats4_4.1.0                SingleCellExperiment_1.16.0 ROCR_1.0-11                 robustbase_0.99-0           tensor_1.5                 
 [26] listenv_0.9.0               MatrixGenerics_1.6.0        labeling_0.4.3              GenomeInfoDbData_1.2.7      numbers_0.8-5              
 [31] polyclip_1.10-4             farver_2.1.1                rprojroot_2.0.3             parallelly_1.36.0           vctrs_0.6.3                
 [36] generics_0.1.3              xfun_0.40                   timechange_0.2.0            diptest_0.76-0              R6_2.5.1                   
 [41] doParallel_1.0.17           GenomeInfoDb_1.30.1         clue_0.3-64                 flexmix_2.3-19              bitops_1.0-7               
 [46] spatstat.utils_3.0-3        DelayedArray_0.20.0         promises_1.2.1              scales_1.2.1                nnet_7.3-19                
 [51] gtable_0.3.4                globals_0.16.2              goftest_1.2-3               rlang_1.1.1                 GlobalOptions_0.1.2        
 [56] splines_4.1.0               lazyeval_0.2.2              spatstat.geom_3.2-5         yaml_2.3.7                  reshape2_1.4.4             
 [61] abind_1.4-5                 httpuv_1.6.11               tools_4.1.0                 ellipsis_0.3.2              RColorBrewer_1.1-3         
 [66] BiocGenerics_0.40.0         ggridges_0.5.4              Rcpp_1.0.11                 plyr_1.8.8                  zlibbioc_1.40.0            
 [71] RCurl_1.98-1.12             deldir_1.0-9                pbapply_1.7-2               GetoptLong_1.0.5            cowplot_1.1.1              
 [76] S4Vectors_0.32.4            zoo_1.8-12                  SummarizedExperiment_1.24.0 ggrepel_0.9.3               cluster_2.1.4              
 [81] here_1.0.1                  magrittr_2.0.3              magick_2.7.4                data.table_1.14.8           scattermore_1.2            
 [86] circlize_0.4.15             lmtest_0.9-40               RANN_2.6.1                  fitdistrplus_1.1-11         matrixStats_1.0.0          
 [91] hms_1.1.3                   patchwork_1.1.3             mime_0.12                   evaluate_0.21               xtable_1.8-4               
 [96] mclust_6.0.0                IRanges_2.28.0              gridExtra_2.3               shape_1.4.6                 compiler_4.1.0             
[101] KernSmooth_2.23-22          crayon_1.5.2                htmltools_0.5.6             later_1.3.1                 tzdb_0.4.0                 
[106] ComplexHeatmap_2.10.0       MASS_7.3-60                 fpc_2.2-10                  Matrix_1.5-1                cli_3.6.1                  
[111] parallel_4.1.0              igraph_1.5.1                GenomicRanges_1.46.1        pkgconfig_2.0.3             sp_2.0-0                   
[116] plotly_4.10.2               spatstat.sparse_3.0-2       foreach_1.5.2               XVector_0.34.0              digest_0.6.33              
[121] sctransform_0.3.5           RcppAnnoy_0.0.21            spatstat.data_3.0-1         rmarkdown_2.24              leiden_0.4.3               
[126] uwot_0.1.16                 kernlab_0.9-32              shiny_1.7.5                 modeltools_0.2-23           rjson_0.2.21               
[131] lifecycle_1.0.3             nlme_3.1-163                jsonlite_1.8.7              BiocNeighbors_1.12.0        limma_3.50.3               
[136] viridisLite_0.4.2           fansi_1.0.4                 pillar_1.9.0                lattice_0.21-8              fastmap_1.1.1              
[141] httr_1.4.7                  DEoptimR_1.1-2              survival_3.5-7              glue_1.6.2                  FNN_1.1.3.2                
[146] png_0.1-8                   prabclus_2.3-2              iterators_1.0.14            bluster_1.4.0               class_7.3-22               
[151] stringi_1.7.12              irlba_2.3.5.1               future.apply_1.11.0